Frequency-dependent fitness induces multistability in coevolutionary dynamics 
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Evolution is simultaneously driven by a number of processes such as mutation, competition and 
random sampling. Understanding which of these processes is dominating the collective evolution- 
ary dynamics in dependence on system properties is a fundamental aim of theoretical research. 
Recent works quantitatively studied coevolutionary dynamics of competing species with a focus 
on linearly frequency-dependent interactions, derived from a game-theoretic viewpoint. However, 
several aspects of evolutionary dynamics, e.g. limited resources, may induce effectively nonlinear 
frequency dependencies. Here we study the impact of nonlinear frequency dependence on evolu- 
tionary dynamics in a model class that covers linear frequency dependence as a special case. We 
focus on the simplest non-trivial setting of two genotypes and analyze the co-action of nonlinear 
frequency dependence with asymmetric mutation rates. We find that their co-action may induce 
novel metastable states as well as stochastic switching dynamics between them. Our results reveal 
how the different mechanisms of mutation, selection and genetic drift contribute to the dynamics 
1 and the emergence of metastable states, suggesting that multistability is a generic feature in systems 

with frequency-dependent fitness. 
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Selection, random genetic drift and mutations are the processes underlying Darwinian evolution. For a long time 
population geneticists have analyzed the dynamics in the simplest setting consisting of two genotypes evolving under 
these processes In those studies, a genotype represents an individual's genetic makeup, completely determining 
all relevant properties of the individual. A key concept is the so-called fitness of a genotype which represents the 
' selection pressure for the individuals. The fitness defines the expected number of offspring an individual will produce. 
Thus, selection acts on fitness differences preferring individuals with higher fitness over individuals with lower fitness. 
Usually it is assumed that individuals have fixed fitnesses defined by their genotype alone fij-fcj. Yet, experimental 
studies have revealed that many natural systems exhibit frequency-dependent selection ks46|, which means that an 

■ individual's fitness not only depends on its genotype, but also on its interactions with other individuals and hence on 
the frequency of the different genotypes in the population. Although such frequency-dependent selection had already 

■ been studied early by Crow and Kimura 0, only recently has it received more attention 043 ■ I n these theoretical 
and computational studies, individuals' interactions are represented by interaction matrices from game theory. This 

. leads to a frequency dependence where the fitness depends directly on the interaction parameters in a linear way. 
However, fitness may depend on many diverse factors such as cooperation (i.e. individuals acting together to increase 
their fitness [12, EH) and resource competition, so that certain systems may exhibit frequency-dependent fitness 
that is nonlinear. For example, in experiments certain hermaphrodites exhibit such nonlinear fitness-dependence @. 
To the best of our knowledge the impact of such nonlinear dependencies on coevolutionary dynamics has not been 
investigated theoretically. 

In this article we show that nonlinear frequency dependence @, [l5| may induce new stable configurations of the 
evolutionary dynamics. Furthermore, we study the impact of asymmetric mutation probabilities on the dynamics 
[3 E3> which was also neglected in most models until now 0, M, O EH- As in previous works on coevolutionary 
dynamics we base our work on the Moran process in a non-spatial environment which is a well established model to 
study evolutionary dynamics and was already used in many applications [l8| . The Moran process is a stochastic birth- 
death process which keeps the population size constant [19j. Therefore, in a two-genotype model the system becomes 
effectively one-dimensional, so that the dynamics may be described by a one-dimensional Markov chain with transition 
rates defined by the Moran process. We derive the stationary probability distribution of the system dynamics via 
its Fokker-Planck equation [20j. Sharp maxima of the distribution reveal metastable points of the dynamics and a 
multitude of such maxima lead to stochastic switching dynamics between multiple stable points. 

The article is structured as follows. In Section II we introduce the model details and in Section III we derive the 
Fokker-Planck equation describing the probabilistic dynamics of the population. Using this equation we derive the 
stationary probability distribution that describes the long-time behavior of the system. In Section IV we analyze this 
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probability distribution, which yields information about the impact of nonlinear frequency-dependent selection and 
of different mutation probabilities on the coevolutionary dynamics. In Section V we give a summary and discuss our 
results. 

II. THE MODEL 

Consider a population of N individuals evolving in a non-spatial environment, meaning that the population is 
well-mixed so that each individual interacts with all other individuals at all times. In this population the individuals 
may assume one out of two genotypes A and B. The population sizes &a and ks (kA + fcs = N) evolve according 
to the time-continuous Moran process described in the following, cf. fl9j |. The number of individuals Ua of genotype 
A completely determines the state of the system as kg = N — k A - At all times the interactions of the individuals 
determine the actual (frequency-dependent) fitness, so that an individual's fitness of genotype A or B is defined by 
a fitness function /a(^a) or /b(&a) respectively. The fitness functions /a(^a) and /b^a) may be any functions 
with the only condition that fi(kA) > for all fc^ G [0,N] as negative fitness is not defined. At rate /a(&a)&A an 
individual of type A produces an identical offspring which may mutate to genotype B with probability hab- This 
applies analogously to genotype B. Then one individual of the population is chosen randomly to die, so that the 
population size N stays constant and the variables ki change maximally by 1. This is the so-called Moran process 
which was originally introduced for a population of two genotypes with fixed fitnesses and no mutations occurring, 
cf. [llj]. However, this process is easily generalizable to more genotypes and frequency-dependent selection in the way 
described above, cf. [Til [2l|. 

Note that in our model the rate of reproduction of e.g. type A is directly determined by the term /ca/a^a) as e.g. 
inj22|. In other models the fitness is first normalized so that the rate of reproduction is given by kAfA(kA)/f(kA) 
[13j . where 

7(k A ) = ^[k A fA(k A ) + (N- k A )f B {k A )] (1) 

is the population's mean fitness. While in both of these models the events occur with the same probability, the times 
between the events differ by a common factor determined by the mean fitness ([1}. Thus, these models exhibit a 
quantitatively different time course, but the same event sequences. 

Until now, usually linear functions have been considered for fi(kA)- However, in many applications nonlinear 
functions seem more appropriate [(I Il5|. For exam ple, cooperation effects in game theory induce functions linearly 
increasing in kA, so that /a(^a) = 1 + a ■ kA [ill. Il3|. This would mean that fitness increases infinitely with the 
population size of genotype A. Yet, in any habitat there is only a limited amount of resources available for living. 
Therefore, fitness should decline when the population of one genotype becomes too large as all of the individuals 
will compete for the same resources. We conclude that a fitness function including these two effects has to contain a 
nonlinear factor, e.g. 

f A (k A ) = l + a-k A -b-k A (2) 

where a > and b > 0. Here we choose a quadratic nonlinear factor as a simple example. Naturally, all other 
nonlinear functions could be applicable as well. 

Note that linearly increasing fitness functions such as /a^a) = 1 + a - kA and /b(^a) = 1 + b • &A lead to quadratic 
relations for the mean fitness /(^a) (see equation (TTJ) of the population which becomes 

J{k A ) = l + bkA + {a-b) 1 ^. (3) 

The quadratic factor originates from the fact that the linear fitness function in ([1]) is multiplied by the number of 
individuals leading to a quadratic dependence on fc^, e.g. &a/a(&a) = &a(1 + Q^a) = &A + ak\. This factor enters 
the replicator equation in evolutionary game theory, thus leading to a quadratic fitness dependence in standard game 
theoretic problems [111 . |22j . However, this does not imply nonlinear interactions and thus fitness effects such as 
presented in equation ([2]). Therefore, the fitness functions that can be generated by evolutionary game theory are a 
special case of the functions that occur in the theory presented here. Fitness functions as for example presented in 
equation ([2]) will - as we show in this article - lead to dynamics with more complex stability properties. 

Usually, when analyzing the dynamics of such a model system as described above, it was assumed either that no 
mutations occur at all (py = 0) 0, EH Ell ° r that the mutation pro babilities are equal (fiAB = ^ba) [2 0, H3 • 
However, usually mutation probabilities can be highly diverse (l6l. |17|. This has been considered in some examples 
[|| , but until now the effect of different mutation probabilities in the introduced model system has not been studied 
systematically. Here, we explicitly study the effects caused by asymmetric mutation rates (see e.g. figure [5]). 
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III. ANALYSIS 



In the following we assume the fitness of the individuals to be given by 

f A (k) = l + g A (k) 

and 

f B (k) = l + g B (k) 



(4) 



(5) 



with the interaction functions g A > — 1 and g B > —1- The model system is effectively one-dimensional with the 
variable k = k A = N — kg completely describing the system state [12]. At each event k may change by at most ±1 
depending only on the actual state of the system, so that the dynamics are Markovian. Let the transition k — > k + 1 
occur with the rate r k . Then this rate is determined by the fitnesses and mutation probabilities in the following way: 

• Each individual in the population receives offspring with rate fj(k) with j 6 {A, B}. This means that genotype 
A receives one offspring with rate k ■ f A (k). 

• Such an offspring increases the population A with probability 1 — p> AB , and population B in case it mutates 
with probability p A s- 

• One individuum is chosen uniformly to die belonging to genotype B with probability (N — k)/(N + 1). 

• Equivalently the number of individuals of genotype A can increase if genotype B receives one offspring, which 
mutates to genotype A with probability pba, an d one individuum of B is chosen to die. 



Taken these processes together the transition rate for k — > k + 1 is 



rt = (1 - HAB)kf A {k) 



N 



Hba(N -k)f B (k) 



N -k 



(G) 



N + l ^ y ,J v ' N + l 
Decreasing k to k — 1 is equivalent to increasing the number of individuals of genotype B, so that the transition rate 

k k 



{1 - Hba)(N - k)f B (k) 



VABkf A {k) 



N+l ^ " w N+l 
is obtained analogously. The master equation of the Markov chain is then given by 

dp k (t) 



dt 



= rt^Pk-iit) - (r k + r+)p k {t) + r k+lPk+1 (t) k e {0, N} 



(7) 



(8) 



where we define r 



+ 



' N+l 



0. As also r n 



' N 



0, the system has reflective boundaries at k = and k = N. If 



fJ'AB = (/j-ba = 0) then k A = N (k A = 0) is an absorbing state of the dynamics. 

It is possible to derive the exact solution to such an equation similarly to the derivations in (23| which we do in 
Appendix B. However, the obtained solutions are not very explicit and their implications are hard to grasp. It is thus 
more useful to advance to the Fokker-Planck equation [20] describing the system in the limit of large N . As previous 
works have shown, this approximation already leads to very good results for moderate population sizes of N < 1000, 
cf. [13. We use the transformation 



^[0,1], 



t 



p{x,s) = p xN {sN)N, gj(x) 



IMj ■ N, 
9j (xN)-N 



(9) 
(10) 



yielding the Fokker-Planck equation. The scaling factors are derived in Appendix A where the derivation of the 
Fokker-Planck equation in the scaling limit N — > oo is shown in detail (Equation (|3Tj)). This leads to 



dp(x, s) 
ds 



d_ 

dx 



(/2(1 - 2x) + [g A (x) - sb (a)] x(l - x) - Aft) p(x, s) 



d_ 

dx 2 



2 r 



x(l — x)p(x, s) 



with the normalization condition 



p(x, s)dx = 1 



(11) 



(12) 
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for all s > 0. Here §a — NgA(Nx) and gs = NgB(Nx) are the rescaled interaction functions defined in (ITUl) . 
ft = N{[iab + Mba)/2 is the rescaled mean mutation rate and Afi = N([1ab ~ (J-ba)/^ is the rescaled mutation rate 
difference. These functions have to be rescaled to obtain a non-degenerate limit, the so-called weak selection regime. 
Notice, that the drift term 



D (1) {x) = fi(l - 2x) + [g A (x) - g B (x)} x(l - x) - Afi 



(13) 



contains three different effects. The fitness difference at each point x causes a selective drift, the mean mutation rate 
causes a drift directed to x = 1/2 and the mutation rate difference causes a one-directional drift towards one of the 
boundaries x = or x = 1. The diffusion term 



£> (2) (x) =x(l-x) 



(14) 



reflects the undirected genetic drift. 

If fitness differences are on a scale of O (iV _1 ), then mutation, selection and genetic drift all act on the same scale 
[l3L [2^| . Then the dynamics are characterized by an interplay of these different effects. On the other hand, in the 
strong selection regime - where fitness differences are on a scale 0(1) - genetic drift becomes negligible and the 
dynamics in the bulk of the system become deterministic for large N [25fl . However, as the factor x{l — x) vanishes 
for x — ¥ and x — > 1 the effects of mutations play an important role. Clearly, without mutations x = and x — 1 are 
absorbing states, but even in the strong selection regime they become non-absorbing for any positive mutation rate. 
This is reflected in the form of the drift term D^\x), as only the factors reflecting mutations fi(l — 2x) — Afi remain 
non-zero in the limits x — > and x — >■ 1. 

For a one-dimensional Fokker-Planck equation pip the stationary distribution is (20| 



p*(x) =Ce~* {x) 



with the potential 



dx 



DW(x) 

which is defined up to a constant irrelevant to our calculations as it is canceled by the normalization 

1 



C 



Thus, we obtain 



$(x) = In (x(l - x)) - 



[7(1 - 2x) 



(15) 



(16) 



(17) 



x(l 



+ 9a(x) -<}b(x) ~ 



Afi 



dx 



In (x(l — x)) — fi ln(x(l — x)) + Afi In 



x(l — x) 
- I \9a{x) -g B {x)]dx 



where we used 



(1 - 2x) -- 

The stationary solution becomes 

P *{x) = C{x{l~x)t~ 1 



dx 



[x(l-x)\. 



-A/i 



of[9A(x)-g B (x)]dx 



(18) 



(19) 



where the constant C defined in (|17|) has to be calculated numerically for given parameters. 

This stationary distribution contains contributions from the above mentioned four different effects: 

1. Genetic drift is reflected by (x(l — a;))" 1 which diverges for x G {0, 1} at the boundaries. 

2. The mean mutation rate fi > causes the balancing term (x(l — x)Y . 

3. The asymmetry in the mutation probabilities causes the term (x/(l — x))^ A ' 1 which diverges at one boundary 
and vanishes at the other. 

4. All frequency-dependent selection effects are contained in the exponential factor e^ A ^~^ B ^ x '^ dx which can take 
various shapes. 
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Figure 1: The stationary distribution for the model system exhibits 3 maxima corresponding to metastable points, (a) shows 
the theoretical curve from equation (|19|) (red, solid) in perfect agreement with data from simulations with N = 1000 (blue, 
x). (b) shows the fitness functions of genotype A (blue, solid) and B (red, dashed). Interaction parameters are a a = 0.083, 
6a = 0.05, as = 0.177 and fes = 0.2 (see equation Q and (|20[0 while the mutation rate is fi = 0.5 (Ap, = 0). 



IV. MULTIPLE STABLE POINTS 



What are the possible shapes for the stationary distribution p*(x) in the form given by equation (fTTJl) ? Until now, 
the term (x/ (1 — x)) _A ' i representing the asymmetric mutation probabilities has not been considered to our knowledge 
and the interaction functions g A (x) and </_b(x) in the selection term were considered to be at most linear in x [l3j. 
We should therefore be interested in the effects of nonlinear interaction functions and asymmetric mutation rates. 

Let us first analyze the dynamics for nonlinear interaction functions [fj, |l5j describing the effects of cooperation and 
limited resources as described by equation ©, so that 

g A {x) = N (a A x - b A x 2 ) g B (x) = N (a B {l - x) - b B {\ - xf) . (20) 

Already with A/i = such interaction functions induce dynamics stochastically switching between three metastable 
points. An example is shown in Figure[TJ where the theoretically calculated stationary distribution from equation (fT9|) 
is shown together with data from simulations with a population of N = 1000 individuals which is enough to obtain 
almost perfect fitting (cf. also Figure 2]). The fitness functions (Figure [TJd) both show at first an increase on increasing 
the number of individuals of genotype A or B from due to cooperation effects and then a strong decrease due to 
resource competion. As the resulting fitness functions are asymmetric, also the stationary distribution is asymmetric. 
There is a maximum at x — due to genetic drift. Selection drives the dynamics towards a metastable state at 
x w 0.2, because for x < 0.2 genotype A is fitter than B thus increasing in frequency and for x > 0.2 genotype A 
is less fit than B thus decreasing in frequency (cf . Figure Q~|d) . The maximum of the stationary distribution is thus 
exactly at the point where f A {x) = f B (x). For x > 0.7 again genotype A is fitter than B and thus the dynamics 
are driven towards x = l by selection as well as genetic drift. The mutational force induced by /t = 0.5 increases the 
height of the maximum at x « 0.2 driving the system away from the maxima at x = and x = 1. So, in this example 
genetic drift, mutation and selection all significantly influence the dynamics. 

Let us now study the influence of the mutation rates in detail. Interestingly, for asymmetric mutation rates (A/I ^ 0) 
the factor (x/(l — x))~ A '- 1 always diverges in the interval [0, 1]. For A/t > it diverges for x — > 0, otherwise for x — > 1. 
This can cause the emergence of a maximum of the stationary distribution at x = (or resp. x = 1). Figure [5] 
shows an example where due to asymmetric mutation rates a maximum occurs at x = 1, when for A/2 = there is 
an absolute minimum at x — 1 as this stable state minimizes the population's mean fitness. This means, that the 
system dynamics are mutation dominated near x — 1. Furthermore, the maximum located at x = for Afl = is 
shifted for A/t > to values x > causing a minimum of the stationary distribution at x = 0. Thus, in this example 
both selection and the asymmetry in mutation probabilities are the driving forces of the system's dynamics. 

For low mutation rates the population has a high tendency to become dominated by one of the two genotypes for 
long times which is illustrated in Figure Thus, the dynamics stay close to the edges x — or x — 1 waiting for one 
of the few mutations to occur. For low overall mutation rates differences in the mutation rates have only weak effects. 
High mutation rates cause the population to be drawn towards a mixture of both genotypes (x ~ 0.5). Even more, for 
high mutation rates the differences in the mutation rates A/i can have an important influence on the dynamics. For 
instance Figure \3]p illustrates a shift of an existing stable point, which completely vanishes for high A/t. We conclude 
that asymmetric mutation rates can cause the emergence of new stable states (as in Figure [5]) , the disappearance of 
existing ones (as in Figure [3]), and shift existing stable states to new positions. 
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Figure 2: Asymmetric mutation rates cause the emergence of a new maximum of the stationary distribution, (a) shows the 
stationary distribution of the Fokker-Planck equation for Aft — (red, solid) and Aft — —0.5 (gray, dashed) for a system with 
selective advantage for genotype B as shown in (b) (blue, solid: genotype A; red, dashed: genotype 73). Solid and dashed 
lines in (a) show the theoretical curves from equation (|19|) , crosses the data from simulations with TV = 1000 (x: Aft — 0, 
+: Aft = —0.5). Interaction parameters are oa = 0.077, 6a = 0.1, as = 0.2 and 6b = 0.15 (see equation (|20I0 while mean 
mutation rate is ft = 1. 




Figure 3: The mutation rates strongly influence the system's dynamics, (a) shows stationary solutions of an example system 
with interaction functions as in equation (|20[) for three different symmetric (Aft = 0) mutation rates: ft = 1 (blue, solid), 
ft = 0.1 (red, long dashed) and ft — 0.01 (green, short dashed). This illustrates that for low mutation rates the dynamics stay 
close to the system's edges for long times. Here, differences in the mutation rates only slightly affect the shape of the stationary 
distribution. On the other hand, (b) shows that for high mutation rates (ft = 5) the dynamics are drawn stronger towards the 
middle. Here, differences in the mutation rates have a stronger impact, as the three curves with Afl = (blue, solid), Afl — 2.5 
(red, long dashed) and Aft — 4.5 (green, short dashed) demonstrate. The increasing Aft shifts the existing stable state towards 
the edge of the system until it vanishes. System parameters were N — 1000 and a symmetric interaction function according to 
equation (|20[) with one stable state at x = 0.5 with <ja = as = —0.01, 6a = 6_g = 0.005. 



As we saw above, the fit of the theoretically calculated stationary distribution to simulation data is almost perfect 
for weak selection, i.e. the maximal fitness difference satisfies A/ max <C 1, cf. e.g. [13]. To quantify the quality of 
the fit we measured the stationary distribution with the same parameters as above for different population sizes N. 
This is shown in Figure 2] demonstrating that for increasing N the data gets ever closer to the theoretically obtained 
distribution. For this we define the empirical distribution 

« T mcas 

E^+^-fc) (21) 

^meas . q 

where (X t : t > 0) is the discrete process defined in equation ([5]). T m j x is a time large enough for the process to 
reach stationarity and depends on the system size, as does the measurement time T moas which is specified in Figure 
|4l Using this definition we quantify the quality of the fit using the distance measure 

N-l 



k/N+l/2N 

:r A / p*(x)dx 
k/N-l/2N 



(22) 
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Figure 4: The theoretically obtained stationary distribution well fits data from simulations for large population sizes N. (a) 
shows the theoretical curve from equation (|f 9[) (black, solid) for the same system as in Figure [T] Data points show empirical 
distributions as defined in equation (|21[) obtained in simulations for N — 50 (blue, x), N — 100 (red, *), N = 500 (green, •) 
and N = 1000 (gray, +). (b) shows the distances d and d max of simulation data and theoretical curve for different N for the 
mean distance measure d defined in equation (|22p (blue, x ) and the maximum distance measure d max defined in equation (|23[) 
(red, *), demonstrating that the distance decays for increasing N. The solid line N~ 3 ^ 2 and the dashed line 2N~ 1 ^ 2 are added 
as a guide to the eye for the relation between measured distances and system size N. The measured empirical distributions for 
both (a) and (b) were obtained by simulating the system dynamics from an initial state drawn from p*(x) for a mixing time 
Tmix = lOOiV and then recording the density for a time r meas = 10iV 2 . 



which takes the mean distance of the measured empirical distribution irk from the theoretical distribution for every 
point except k = and k = N where the density p*(x) diverges. Here, the the theoretical distribution is calculated by 
integrating the theoretical density over bins of the size 1/N around the points k/N. Figure 0}) shows that the mean 
distance d decays with increasing N. 

The decay of this quantity is dominated by the slow convergence at the domain boundaries k = and k = N. 
Therefore, we additionally defined the maximum distance measure 



dmax '■— max 

fc£[l,JV-l] 



k/N+l/2N 

ir k - I p*(x)dx 



(23) 



which gives the maximum distance between the measured probability distribution from the theoretical distribution 
for all points except k = and k = N due to the same argument as above. Figure |3b shows that the maximum 
distance d max decreases with increasing N, however slower than the mean distance d due to the divergence of p(x) 
at the boundaries. Altogether we conclude, that the theoretical curve fits the data very well already for N £3 1000. 
However, near the domain boundaries special care has to be taken, as divergences of the theoretical curve lead to 
larger deviations. 

Furthermore, as the Fokker-Planck approximation only holds for weak selection we study the quality of the solution 
obtained throught the Fokker-Planck equation in dependence of the selection strength. For this we introduce a scaling 
factor £ to the interaction functions, so that they become 

g A (x) = £N (a A x - b A x 2 ) g B (x) = £N (a B (l - x) - b B {l - xf) . (24) 

We find that the stationary solution of the Fokker-Planck equation well approximates the solution of the Master 
equation 

u k-l r+ 
llj=0 r - 

P* k = (25) 

1^1=0 1 1.3=0 r - 

3 + 1 

for selection strengths up to £ of the order of 1 (cf. Figure [5]). We derive the solution (|25|) of the Master equation in 
Appendix B. Note that for the derivation of this solution no approximation is necessary and it hence fits data from 
simulations perfectly up to an error due to finite sampling in simulations. As Figure [5] illustrates, for strong selection 
£ £3 1 the solution of the Fokker-Planck equation does not fit the exact solution of the Master equation perfectly, but 
still catches the overall trend of the dynamics. The error, quantified by the mean distance measure analogous to (|2"2"1) . 
increases with £ as a power law both in comparison with the solution from the Master equation and with data from 
simulations, cf. Figure [SJi. All in all, for weak selection the Fokker-Planck approximation works well while for strong 




0.1 0.5 1 > 5 10 50 



Figure 5: The stationary distribution ()19|) obtained from the Fokker-Planck equation well fits the exact solution Q25p from the 
Master equation for weak selection, but not for strong selection, (a) shows the stationary solution from the Fokker-Planck 
equation (red, solid) for the system with interaction functions defined by equation (|24[) for very weak selection £ = 0.1 together 
with the solution from the Master equation (blue, x). The distributions of the same system for weak selection £ = 1 and strong 
selection £ = 50 are shown in (b) and (c), respectively. All curves are plotted logarithmically to better allow a comparison of 
the deviations over all scales, (c) shows that the curve from the Fokker-Planck equation does not fit the exact Master solution 
very well. This is quantified in (d) showing the mean distance d between Fokker-Planck and Master solution (blue, x) and 
between Fokker-Planck solution and an empirical density obtained from simulations (red, *). The system parameters were 
N — 1000, ft — 0.5, a a — clb = —0.01 and 6a = &s = 0.005. The measured empirical distributions were obtained by simulating 
the system dynamics from an initial state drawn from p* (x) for a mixing time T m i x = 10 5 and then recording the density for a 
time Tmeas = 10 7 . 



selection it does not perfectly predict the stationary distribution, but still reflects the overall trend of the dynamics. 
If this approximation is not satisfying, the direct solution (|25|) of the Master equation yields the exact distribution 
fitting the data for all selection strengths. 

We have shown above that the dynamics of a two-genotype system can exhibit multiple stable points induced by 
nonlinear selection. Even more, there is no theoretical limit to the number of stable points in the system. Assume 
for example fl — 1 and A/2 = 0, so that mutation exactly balances genetic drift. Then, the stationary distribution 
(|19p has a maximum at each point, where exp \J cia{x) — gB{x)dx\ has a maximum. Theoretically (]a(x) and cib{x) 
can be any function with an arbitrary amount of extreme points in [0,1], so that there is no limit to the stationary 
distribution's number of maxima, if selection dominates the dynamics. However, for finite N the number of possible 
maxima is naturally limited by N/2. Figure [S] shows an example, where we used periodic interaction functions 

<ja{x) = a (1 + sin(/3x)) and c)b(x) = a (1 + cos(/3.t)) . (26) 

Although this is not a realistic interaction function in most applications, it demonstrates what is theoretically possible 
in the introduced system. 



V. SUMMARY 



In this article we analyzed a two-genotype system in a very general setting with (possibly) asymmetric mutation 
probabilities and nonlinear fitness functions jg, in finite populations. The underlying Moran process is a well 
established model P, UH, [l3j to gain an understanding of the interplay of selection, mutation and genetic drift in 
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X t 

Figure 6: Multiple stability dynamics for periodic interaction functions, (a) shows the theoretically computed stationary 
distribution from equation ()19|) (red, solid) together with simulation data with TV = 1000 (blue, x ) for the interaction functions 
given in (|26[) . The computed stationary distribution diverges for x — ¥ and x — ► 1 while the simulation data remains finite 
due to the finite number of individuals, (b) shows a sample path which exhibits switching between the different maxima of the 
distribution. Parameters are a = 30, /3 = 25, ft = 0.8 and Aft = 0. 

evolutionary dynamics. However, the Moran process is studied mostly with symmetric mutation probabilities and at 
most linear interaction functions. We reasoned that neither need mutation probabilities be symmetric - as experiments 
have shown, that mutation probabilities are often asymmetric [l6| - nor can all interaction effects be described by 
linear interaction functions. For example cooperation in game theory leads to an interaction function increasing 
linearly in the frequency of the cooperating genotype. Yet, in many applications also cooperators in the end compete 
for the same type of resource which is limited. Therefore, due to limited resources a population being too large cannot 
be sustained leading to a decrease in the fitness. There is no linear function that can reflect both of these effects at 
the same time. 

We derived the Fokker-Planck equation describing the dynamics of the number of individuals k of genotype A in 
the limit of large population sizes N. We quantified the quality of the Fokker-Planck approach for an example (cf. 
Figure 4) where the difference of simulation data and theoretical solution became almost not detectable for population 
sizes larger than N £3 1000. Actually, if the system exhibits absorbing states then the Fokker-Planck method does 
not work to study the corresponding quasi-stationary distributions. Instead, WKB methods are more appropriate to 
describe the system dynamics as for example in [Icj . where fixation resulting from large fluctuations was studied. In 
our model system no such absorbing states exist, as long as the mutation rates are positive (/i.y 7^ 0) and therefore 
the Fokker-Planck equation is appropriate to describe the system dynamics. 

We identified the individual effects of selection, mean mutation rate and mutation difference as well as genetic 
drift and derived the stationary probability distribution as determined by the Fokker-Planck equation. Analyzing the 
distribution, we found that asymmetries in the mutation probabilities may not only induce the shifting of existing 
stable points of the dynamics to new positions, but also lead to the emergence of new stable points. Thus, a genotype 
that has a selective disadvantage can anyway have a stable dynamical state where its individuals dominate the 
population due to a higher mutational stability (see Figure [2]). Further we found, that dynamic fitness leads to 
multiple stable points of the dynamics induced by selection and also genetic drift. We showed an example (Figure [T]) 
where three stable points exist, two caused by selection and one by genetic drift. 

We conclude that frequency-dependent fitness together with asymmetric mutation rates induces complex evolution- 
ary dynamics, in particular if the interactions imply nonlinear fitness functions. Theoretically, there is no limit for 
the number of stable points that the dynamics can exhibit (see Figure [5]) . All in all, we interprete our results such 
that in real biological systems multiple metastable equilibria may exist, whenever species interact in a way complex 
enough to imply a fitness that nonlinearly depends on frequency. As a consequence, one species may exhibit a certain 
frequency for a long time before a sudden shift occurs and then a new frequency prevails. Such a change may thus 
occur even in the absence of changes of the environment; it may be induced as well by a stochastic switching from 
one metastable state to another due to complex inter-species interactions. 

The Moran process is a standard tool to gain theoretical insights into experimental data [18( . Of course, we do not 
propose here that any experimental setup may be exactly described by the Moran process. However, we think that it 
should be feasible to develop an experiment where two different mutants evolve with asymmetric mutation rates. To 
find an experimental setup where the two genotypes also exhibit nonlinear fitness could however prove more difficult. 
Rather, our study is a theoretical study indicating that nonlinear fitness may be the cause of multible stable states 
when observed in experimental data. 
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For further studies on systems with more genotypes it should be useful to combine our considerations presented 
here with the work of Traulsen et al. (l3| , where an analysis of systems with more than two genotypes was carried 
out. Extending those results it may be possible to gain a better understanding of the effects of nonlinear interactions 
for many different genotypes. Also, it may be interesting to study the effects of changing interactions, where the 
interactions change according to the system dynamics [26]. Thus, our study might serve as a promising starting point 
to investigate how nonlinear frequency dependencies impact evolutionary dynamics in complex environments. 
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VI. APPENDIX A THE FOKKER-PLANCK EQUATION 

The master equation ([S} may be transformed to a Fokker-Planck equation in the limit of large N j2^|. We introduce 
the transformation 



x = - s = t- F(N) 



jlij = ■ G{N) 



together with the rescaled functions 



p(x, s) = p Nx {t)N g 3 {x) = g 3 {k) ■ H(N). 



(27) 



(28) 



We fix the scaling functions F(N), G(N) and H(N) such that in the limit N — > oo all terms in equation ([5]) remain 
finite so that mutation, selection and genetic drift all act on the same scale. We further define 



1 , 1 

— x+ = 1 H 

N + N 



(29) 



and substituting all this into the master equation ([SJ, we obtain 



dp{x, s) 
ds 



F(N) 



N 2 



N 



- {[(1 - //ab)(1 + 9a{x-))x-(1 - x-) + pba{1 + g B {x-)){l - x_) 2 ] p(x-,s) 



[(1 - pba)(1 + 9b(x+)){1 - x+)x + + habO- +9a(x+))x+] p{x+,s) 
- [(1 - /Uab)(1 + 9a(x))x(1 -x) + pba{^ + 9b{x))(1 - x) 2 
+ (1 - Mba)(1 + 9b{x))x{1 -x) + pab(^ + 9a(x))x 2 ] p(x, s)} 



(30) 



We choose F(N) = 1/(N + 1) and G(N) = H(N) — N so that in the limit N -> oo the terms stay finite. Further we 
introduce the mean mutation rate p = N(pab + Msa)/2 and the mutation rate difference A/2 = N(pab — ^ba)/ 1 ^- 
To not overload the notation we drop the time argument s of p(x, s) in the following calculation. This leads to 



dp(x) 
ds 



N 2 {-2x(l - x)p(x) + x+(l - x + )p(x + ) + x-(l - x-)p{xJ] 



N 



-(1 - 2x) 2 p(x) + i(l - 2x+) 2 p{x+) + ~(1 - 2x-fp{xJ) 



+N {g A (x^)x-(l -xJ)p{x-) ~g A {x)x(l - x)p(x) + g B (x+)x + (l -x + )p(x + ) - g B {x)x(l - x)p(x) 

+ | [(1 - 2x-)p{x^) - (1 - 2x+)p{x+)] 

+Ajl [x + p(x + ) — xp{x) — (1 — X-)p(x-) + (1 - x)p(x)] 

+ [(9a(x+)x 2 + +g B {x+)(x 2 + - x + )) p(x + ) - (g A {x)x 2 +g B {x){x 2 - x)) p(x) 

+ (9a(x-)(x 2 _ - x-) +g B {xJ)(l - x-) 2 ) p(x-) - (g A (x)(x 2 ~ x) + g B (x)(l - x) 2 ) p{x)] 

+ [(9a{x+)x 2 + - g B (x+)(x 2 + - x + )) p{x+) - (gA(x)x 2 - ()b{x)(x 2 - x)) p{x) 

+ (9a(x~)(x 2 _ - x_) - g B (x-)(l - x^) 2 ) p(x_) - (g A {x)(x 2 - x) - g B (x)(l - x) 2 ) p(x)] } 

In the limit N — > oo the different terms with N 2 in front become second order derivatives with respect to x, while 
the other terms become first order derivatives. The terms which have a l/N factor vanish in the limit N — > oo and 
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thus the above equation becomes 

dp d cP 1 

■£ = ~fa W - 2x )p{ x ) + (IM^O - 9b{x)] x(1 - x)p(x)) - Afip{x)} + [x(l - x)p(x)} (31) 
which is the Focker-Planck-Equation of the system. 

VII. APPENDIX B - STATIONARY SOLUTION OF THE MASTER EQUATION 

We directly derive the stationary solution "p* k of the master equation ((5J) using the detailed balance equation 

r t-iP*k-i = r kP*k ( 32 ) 

which applies to any chain with only nearest neighbour transitions (20| . cf. also [23j . Thus, using the rewritten balance 
equation 

Pt = J ^-Pt-i (33) 
r k 

iteratively, we obtain 

Pt=Pof[^-- ( 34 ) 

Finally, we may use the normalization condition 

£p? = 1 (35) 
;=o 

of the stationary distribution to eliminate the factor pg . We then obtain the exact stationary solution of the master 
equation 

_T+_ 

llj=0 

Pt = 2± V" (36) 

1^1=0 llj=0 T - 

which can be evaluated numerically for any transition rates and . For more details on the exact solution of the 
Master equation in the Moran process see for example the work by Claussen and Traulsen (23| . 



3=0 T 3 + i 



N 
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